######################################################
# CreatePopulationVariables.R
# This script creates school-specific population density estimates 
# using WorldPop data. Note that the code for the raster merge
# cannot be executed without coordinates (which have been removed to comply with IRB)
# The raster merge code has been included below for reference however. 
# note this file requires the working directory to be set to "./Merge and Clean Data/"
# Contact Ryan Jablonski, r.s.jablonski@lse.ac.uk with questions
#
# Log
# Created 2018
# Edited for APSR replication 17 August 2023 by Ryan Jablonski
######################################################




rm(list = ls())


rm(list=ls(all=TRUE))
library(groundhog)

groundhog.library(maps,'2023-04-01')

library(foreign)
groundhog.library(sp,'2023-04-01')
library(maptools)
gpclibPermit()
library(classInt)
groundhog.library(rgeos,'2023-04-01')
library(dismo)
groundhog.library(rgdal,'2023-04-01')
library(raster)


rm(list=ls(all=TRUE))
# 
# p4s = CRS("+proj=longlat +datum=WGS84")
# r.proj=CRS("+proj=utm +zone=36 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0")
# 
# #from https://hub.worldpop.org/geodata/summary?id=123
# r=raster(".\\input\\MWI_pph_2015_v2.tif")
# 
schooldata.lc=read.csv("./input/school_treatment_information_lc_all_schools.csv", stringsAsFactors=FALSE)
schooldata.mp=read.csv("./input/school_treatment_information_mp_all_schools.csv", stringsAsFactors=FALSE)
# 
# 
# school.points <- SpatialPointsDataFrame(coords = schooldata.lc[, c("school_longitude","school_latitude")], data = schooldata.lc,
#                                         proj4string = p4s)
# school.points <- spTransform(school.points, r.proj)
# 
# data.extract=extract(r, school.points)
# r.df=(as.data.frame(r))
# 
# data.malawi.merge = cbind(as.data.frame(school.points), data.frame(pop_per_hectacre=data.extract))
# 
# 
# write.csv(data.malawi.merge[,c("school_id", "pop_per_hectacre")], "./input/worldpop_school_lc.csv")

pop.lc=read.csv("./input/worldpop_school_lc.csv")
schooldata.lc$pop_per_hectacre=pop.lc[match(schooldata.lc$school_id, pop.lc$school_id),"pop_per_hectacre"]

write.csv(schooldata.lc, "./output/Schools.forLC.withpop.csv")


# #MP
# 
# school.points <- SpatialPointsDataFrame(coords = schooldata.mp[, c("school_longitude","school_latitude")], data = schooldata.mp,
#                                         proj4string = p4s)
# school.points <- spTransform(school.points, r.proj)
# 
# data.extract=extract(r, school.points)
# r.df=(as.data.frame(r))
# 
# data.malawi.merge = cbind(as.data.frame(school.points), data.frame(pop_per_hectacre=data.extract))
# 
# write.csv(data.malawi.merge[,c("school_id", "pop_per_hectacre")], "./input/worldpop_school_mp.csv")

pop.mp=read.csv("./input/worldpop_school_mp.csv")
schooldata.mp$pop_per_hectacre=pop.mp[match(schooldata.mp$school_id, pop.mp$school_id),"pop_per_hectacre"]

write.csv(schooldata.mp, "./output/Schools.forMP.withpop.csv")


